/* ========================================================================== */
/* === Demo/cholmod_demo ==================================================== */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Demo Module.  Copyright (C) 2005-2013, Timothy A. Davis
 * The CHOLMOD/Demo Module is licensed under Version 2.0 of the GNU
 * General Public License.  See gpl.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * -------------------------------------------------------------------------- */

/* Read in a matrix from a file, and use CHOLMOD to solve Ax=b if A is
 * symmetric, or (AA'+beta*I)x=b otherwise.  The file format is a simple
 * triplet format, compatible with most files in the Matrix Market format.
 * See cholmod_read.c for more details.  The readhb.f program reads a
 * Harwell/Boeing matrix (excluding element-types) and converts it into the
 * form needed by this program.  reade.f reads a matrix in Harwell/Boeing
 * finite-element form.
 *
 * Usage:
 *	cholmod_demo matrixfile
 *	cholmod_demo < matrixfile
 *
 * The matrix is assumed to be positive definite (a supernodal LL' or simplicial
 * LDL' factorization is used).
 *
 * Requires the Core, Cholesky, MatrixOps, and Check Modules.
* Optionally uses the Partition and Supernodal Modules.
 * Does not use the Modify Module.
 *
 * See cholmod_simple.c for a simpler demo program.
 *
 * cholmod_demo is the same as cholmod_l_demo, except for the size of the
 * basic integer (int vs SuiteSparse_long)
 */

#include <ctype.h>
#include <limits.h>
#include <float.h>
#include <signal.h>
#include <stdio.h>

#include "cholmod_demo.h"
#define NTRIALS 100

/* ff1 and ff2 are global variables so that they can be closed by my_handler */
FILE *ff1 ;
FILE *ff2 ;

int my_atof(char *string, double *pnumber)
{
        /* Convert char string to double data type. */
        double          retval;
        double          one_tenth = 0.1;
        double          ten = 10.0;
        double          zero = 0.0;
        int             found_digits = 0;
        int             is_negative = 0;
        char           *num;

        /* Check pointers. */
        if (string == 0) {
                *pnumber = zero;
                return 0;
        }
        retval = zero;

        num = string;

        /* Advance past white space. */
        while (isspace(*num))
                num++;

        /* Check for sign. */
        if (*num == '+')
                num++;
        else if (*num == '-') {
                is_negative = 1;
                num++;
        }
        /* Calculate the integer part. */
        while (isdigit(*num)) {
                found_digits = 1;
                retval *= ten;
                retval += *num - '0';
                num++;
        }

        /* Calculate the fractional part. */
        if (*num == '.') {
                double          scale = one_tenth;
                num++;
                while (isdigit(*num)) {
                        found_digits = 1;
                        retval += scale * (*num - '0');
                        num++;
                        scale *= one_tenth;
                }
        }
        /* If this is not a number, return error condition. */
        if (!found_digits) {
                *pnumber = zero;
                return 0;
        }
        /* If all digits of integer & fractional part are 0, return 0.0 */
        if (retval == zero) {
                *pnumber = zero;
                return 1;               /* 
                                         * Not an error condition, and no
                                         * need to
                                         * continue. */
        }
        /* Process the exponent (if any) */
        if ((*num == 'e') || (*num == 'E')) {
                int             neg_exponent = 0;
                int             get_out = 0;
                long            index;
                long            exponent = 0;
                double          getting_too_big = DBL_MAX * one_tenth;
                double          getting_too_small = DBL_MIN * ten;

                num++;
                if (*num == '+')
                        num++;
                else if (*num == '-') {
                        num++;
                        neg_exponent = 1;
                }
                /* What if the exponent is empty?  Return the current
                 * result. */
                if (!isdigit(*num)) {
                        if (is_negative)
                                retval = -retval;

                        *pnumber = retval;

                        return (1);
                }
                /* Convert char exponent to number <= 2 billion. */
                while (isdigit(*num) && (exponent < LONG_MAX / 10)) {
                        exponent *= 10;
                        exponent += *num - '0';
                        num++;
                }

                /* Compensate for the exponent. */
                if (neg_exponent) {
                        for (index = 1; index <= exponent && !get_out; index++)
                                if (retval < getting_too_small) {
                                        get_out = 1;
                                        retval = DBL_MIN;
                                } else
                                        retval *= one_tenth;
                } else
                        for (index = 1; index <= exponent && !get_out; index++) {
                                if (retval > getting_too_big) {
                                        get_out = 1;
                                        retval = DBL_MAX;
                                } else
                                        retval *= ten;
                        }
        }
        if (is_negative)
                retval = -retval;

        *pnumber = retval;

        return (1);
}

/* halt if an error occurs */
static void my_handler (int status, const char *file, int line,
                        const char *message)
{
        printf ("cholmod error: file: %s line: %d status: %d: %s\n",
                file, line, status, message) ;
        if (status < 0)
        {
                if (ff1 != NULL) fclose (ff1) ;
                exit (0) ;
        }
}

int main (int argc, char **argv)
{
        double resid [4], t, ta, tf, ts [3], tot, bnorm, xnorm, anorm, rnorm, fl,
                anz, axbnorm, rnorm2, resid2, rcond ;
        FILE *f_1, *f_2 ;        
        cholmod_sparse *A , *B_MTR = NULL;
        cholmod_dense *X = NULL, *B, *W, *R ;
        double one [2], zero [2], minusone [2], beta [2], xlnz ;
        cholmod_common Common1, Common2, *cm, *cm1, *cm2 ;
        cholmod_factor *L ;
        double *Bx, *Rx, *Xx ;
        int i, n, isize, xsize, ordering, xtype, s, ss, lnz ;
        int trial, method, L_is_super ;
        int ver [3] ;

        /* ---------------------------------------------------------------------- */
        /* get the file containing the input matrix */
        /* ---------------------------------------------------------------------- */

        ff1 = NULL ;
        ff2 = NULL ;
        if (argc > 1)
        {
                if ((f_1 = fopen (argv[1], "r")) == NULL)
                {
                        my_handler (CHOLMOD_INVALID, __FILE__, __LINE__,
                                    "unable to open matrix A") ;
                }
	
                if ((f_2 = fopen (argv[2], "r")) == NULL)
                {
                        my_handler (CHOLMOD_INVALID, __FILE__, __LINE__,
                                    "unable to open matrix B") ;
                }
        
                ff1 = f_1 ;
                ff2 = f_2 ;
        }
        else
        {
                f_1 = stdin ;
                f_2 = stdin ;
        }
        
        /* ---------------------------------------------------------------------- */
        /* start CHOLMOD and set parameters */
        /* ---------------------------------------------------------------------- */
        cm1 = &Common1 ;
        cholmod_start (cm1) ;
        cm = cm1;
        
        //cm2 = &Common2 ;
        //cholmod_start (cm2) ;
        
        /* use default parameter settings, except for the error handler.  This
         * demo program terminates if an error occurs (out of memory, not positive
         * definite, ...).  It makes the demo program simpler (no need to check
         * CHOLMOD error conditions).  This non-default parameter setting has no
         * effect on performance. */
        cm1->error_handler = my_handler ;
        //cm2->error_handler = my_handler ;

        /* Note that CHOLMOD will do a supernodal LL' or a simplicial LDL' by
         * default, automatically selecting the latter if flop/nnz(L) < 40. */

        /* ---------------------------------------------------------------------- */
        /* create basic scalars */
        /* ---------------------------------------------------------------------- */

        zero [0] = 0 ;
        zero [1] = 0 ;
        one [0] = 1 ;
        one [1] = 0 ;
        minusone [0] = -1 ;
        minusone [1] = 0 ;
        beta [0] = 1e-6 ;
        beta [1] = 0 ;

        /* ---------------------------------------------------------------------- */
        /* read in a matrix */
        /* ---------------------------------------------------------------------- */

        printf ("\n---------------------------------- cholmod_HSTA:\n") ;
        cholmod_version (ver) ;
        printf ("cholmod version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
        SuiteSparse_version (ver) ;
        printf ("SuiteSparse version %d.%d.%d\n", ver [0], ver [1], ver [2]) ;
        A = cholmod_read_sparse (f_1, cm1) ;
        if (ff1 != NULL)
        {
                fclose (ff1) ;
                ff1 = NULL ;
        }
        anorm = cholmod_norm_sparse (A, 0, cm1) ;
        xtype = A->xtype ;
        printf ("xtype = %d\n", xtype) ;
        printf ("norm (A,inf) = %g\n", anorm) ;
        printf ("norm (A,1)   = %g\n", cholmod_norm_sparse (A, 1, cm1)) ;
        cholmod_print_sparse (A, "A", cm1) ;        
        
        if (A->nrow > A->ncol)
        {
                /* Transpose A so that A'A+beta*I will be factorized instead */
                cholmod_sparse *C1 = cholmod_transpose (A, 2, cm1) ;
                cholmod_free_sparse (&A, cm1) ;
                A = C1 ;
                printf ("transposing input matrix\n") ;
        }
        
        FILE* f_;
        if ((f_ = fopen ("Results_", "w+")) == NULL){
                printf("Oh no!\n");
                exit(1);
        } 
        cholmod_write_sparse (f_, A, B_MTR, NULL, cm1) ;
        
        //Bx = (double*)malloc(sizeof(double)*n);
        //keep the following for both cases
        n = A->nrow ; 
        B = cholmod_zeros (n, 1, xtype, cm1) ;
        Bx = B->x ;
        
        //Reading matrix B, which is a dense matrix        
        char* temp;
        double* tmp = (double*)malloc(sizeof(double));
        char line[65532];
        fgets(line, 65532, ff2);
        //printf("%s\n", line);
        temp = strtok(line, " \t");
        //if(temp != NULL)
        //printf("%s\n", temp);
        
        my_atof(temp, tmp);
        Bx[0] = *tmp;
        for(i = 1; i < n; i++){
                temp = strtok(NULL, " \t\n");
                //if(temp != NULL)
                //printf("%s\n", temp);
                my_atof(temp, tmp);
                Bx[i] = *tmp;
        }
        
        FILE* f__;
        if ((f__ = fopen ("Results__", "w+")) == NULL){
                printf("Oh no!\n");
                exit(1);
        }
        cholmod_write_dense (f__, B, NULL, cm1) ;
        
        //Bx = B->x ;
        
        for(i = 0; i < n; i++){
                printf("%.5f\n", ((double*)B->x)[i]);
                printf("%.5f\n", Bx[i]);
        }
        
        
        /* ---------------------------------------------------------------------- */
        /* create an arbitrary right-hand-side */
        /* ---------------------------------------------------------------------- */

/* #if GHS */
/*         { */
/*                 /\* b = A*ones(n,1), used by Gould, Hu, and Scott in their experiments *\/ */
/*                 cholmod_dense *X0 ; */
/*                 X0 = cholmod_ones (A->ncol, 1, xtype, cm) ; */
/*                 cholmod_sdmult (A, 0, one, zero, X0, B, cm) ; */
/*                 cholmod_free_dense (&X0, cm) ; */
/*         } */
/* #else */
/*         if (xtype == CHOLMOD_REAL) */
/*         { */
/*                 /\* real case *\/ */
/*                 for (i = 0 ; i < n ; i++) */
/*                 { */
/*                         double x = n ; */
/*                         Bx [i] = 1 + i / x ; */
/*                         printf("%.5f\n", Bx[i]); */
/*                 } */
/*         } */
/*         else */
/*         { */
/*                 /\* complex case *\/ */
/*                 for (i = 0 ; i < n ; i++) */
/*                 { */
/*                         double x = n ; */
/*                         Bx [2*i  ] = 1 + i / x ;		/\* real part of B(i) *\/ */
/*                         Bx [2*i+1] = (x/2 - i) / (3*x) ;	/\* imag part of B(i) *\/ */
/*                 } */
/*         } */
/* #endif */
    
/*         FILE* f__; */
/*         if ((f__ = fopen ("Results__", "w+")) == NULL){ */
/*                 printf("Oh no!\n"); */
/*                 exit(1); */
/*         } */
/*         cholmod_write_dense (f__, B, NULL, cm1) ; */
        
        cholmod_print_dense (B, "B", cm1) ;
        bnorm = cholmod_norm_dense (B, 0, cm1) ;	/* max norm */
        printf ("bnorm %g\n", bnorm) ;

        /* ---------------------------------------------------------------------- */
        /* analyze and factorize */
        /* ---------------------------------------------------------------------- */

        t = CPUTIME ;
        L = cholmod_analyze (A, cm1) ;
        ta = CPUTIME - t ;
        ta = MAX (ta, 0) ;

        printf ("Analyze: flop %g lnz %g\n", cm1->fl, cm1->lnz) ;

        if (A->stype == 0)
        {
                printf ("Factorizing A*A'+beta*I\n") ;
                t = CPUTIME ;
                cholmod_factorize_p (A, beta, NULL, 0, L, cm1) ;
                tf = CPUTIME - t ;
                tf = MAX (tf, 0) ;
        }
        else
        {
                printf ("Factorizing A\n") ;
                t = CPUTIME ;
                cholmod_factorize (A, L, cm1) ;
                tf = CPUTIME - t ;
                tf = MAX (tf, 0) ;
        }

        cholmod_print_factor (L, "L", cm1) ;

        /* determine the # of integers's and reals's in L.  See cholmod_free */
        if (L->is_super)
        {
                s = L->nsuper + 1 ;
                xsize = L->xsize ;
                ss = L->ssize ;
                isize =
                        n	/* L->Perm */
                        + n	/* L->ColCount, nz in each column of 'pure' L */
                        + s	/* L->pi, column pointers for L->s */
                        + s	/* L->px, column pointers for L->x */
                        + s	/* L->super, starting column index of each supernode */
                        + ss ;	/* L->s, the pattern of the supernodes */
        }
        else
        {
                /* this space can increase if you change parameters to their non-
                 * default values (cm1->final_pack, for example). */
                lnz = L->nzmax ;
                xsize = lnz ;
                isize =
                        n	/* L->Perm */
                        + n	/* L->ColCount, nz in each column of 'pure' L */
                        + n+1	/* L->p, column pointers */
                        + lnz	/* L->i, integer row indices */
                        + n	/* L->nz, nz in each column of L */
                        + n+2	/* L->next, link list */
                        + n+2 ;	/* L->prev, link list */
        }

        /* solve with Bset will change L from simplicial to supernodal */
        rcond = cholmod_rcond (L, cm1) ;
        L_is_super = L->is_super ;

        /* ---------------------------------------------------------------------- */
        /* solve */
        /* ---------------------------------------------------------------------- */

        FILE* f = fopen("Results", "w+");
        
        for (method = 0 ; method <= 0 ; method++)
        {
                double x = n ;

                if (method == 0)
                {
                        /* basic solve, just once */
                        t = CPUTIME ;
                        cholmod_write_sparse (f, A, B_MTR, NULL, cm1) ;
                        cholmod_write_dense (f, B, NULL, cm1) ;
                        X = cholmod_solve (CHOLMOD_A, L, B, cm1) ;
                        cholmod_write_dense (f, X, NULL, cm1) ;
                        ts [0] = CPUTIME - t ;
                        ts [0] = MAX (ts [0], 0) ;
                }
                else if (method == 1)
                {
                        /* basic solve, many times, but keep the last one */
                        t = CPUTIME ;
                        for (trial = 0 ; trial < NTRIALS ; trial++)
                        {
                                cholmod_free_dense (&X, cm1) ;
                                Bx [0] = 1 + trial / x ;        /* tweak B each iteration */
                                X = cholmod_solve (CHOLMOD_A, L, B, cm1) ;
                        }
                        ts [1] = CPUTIME - t ;
                        ts [1] = MAX (ts [1], 0) / NTRIALS ;
                }
                else if (method == 2)
                {
                        /* solve with reused workspace */
                        cholmod_dense *Ywork = NULL, *Ework = NULL ;
                        cholmod_free_dense (&X, cm1) ;

                        t = CPUTIME ;
                        for (trial = 0 ; trial < NTRIALS ; trial++)
                        {
                                Bx [0] = 1 + trial / x ;        /* tweak B each iteration */
                                cholmod_solve2 (CHOLMOD_A, L, B, NULL, &X, NULL,
                                                &Ywork, &Ework, cm1) ;
                        }
                        cholmod_free_dense (&Ywork, cm1) ;
                        cholmod_free_dense (&Ework, cm1) ;
                        ts [2] = CPUTIME - t ;
                        ts [2] = MAX (ts [2], 0) / NTRIALS ;
                }
                else
                {
                        /* solve with reused workspace and sparse Bset */
                        cholmod_dense *Ywork = NULL, *Ework = NULL ;
                        cholmod_dense *X2 = NULL, *B2 = NULL ;
                        cholmod_sparse *Bset, *Xset = NULL ;
                        int *Bsetp, *Bseti, *Xsetp, *Xseti, xlen, j, k, *Lnz ;
                        double *X1x, *X2x, *B2x, err ;
                        FILE *timelog = fopen ("timelog.m", "w") ;
                        if (timelog) fprintf (timelog, "results = [\n") ;

                        B2 = cholmod_zeros (n, 1, xtype, cm1) ;
                        B2x = B2->x ;

                        Bset = cholmod_allocate_sparse (n, 1, 1, FALSE, TRUE, 0,
                                                        CHOLMOD_PATTERN, cm1) ;
                        Bsetp = Bset->p ;
                        Bseti = Bset->i ;
                        Bsetp [0] = 0 ;     /* nnz(B) is 1 (it can be anything) */
                        Bsetp [1] = 1 ;
                        resid [3] = 0 ;

                        for (i = 0 ; i < MIN (100,n) ; i++)
                        {
                                /* B (i) is nonzero, all other entries are ignored
                                   (implied to be zero) */
                                Bseti [0] = i ;
                                if (xtype == CHOLMOD_REAL)
                                {
                                        B2x [i] = 3.1 * i + 0.9 ;
                                }
                                else
                                {
                                        B2x [2*i  ] = i + 0.042 ;
                                        B2x [2*i+1] = i - 92.7 ;
                                }

                                /* first get the entire solution, to compare against */
                                cholmod_solve2 (CHOLMOD_A, L, B2, NULL, &X, NULL,
                                                &Ywork, &Ework, cm1) ;

                                /* now get the sparse solutions; this will change L from
                                   supernodal to simplicial */

                                if (i == 0)
                                {
                                        /* first solve can be slower because it has to allocate
                                           space for X2, Xset, etc, and change L.
                                           So don't time it */
                                        cholmod_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
                                                        &Ywork, &Ework, cm1) ;
                                }

                                t = CPUTIME ;
                                for (trial = 0 ; trial < NTRIALS ; trial++)
                                {
                                        /* solve Ax=b but only to get x(i).
                                           b is all zero except for b(i).
                                           This takes O(xlen) time */
                                        cholmod_solve2 (CHOLMOD_A, L, B2, Bset, &X2, &Xset,
                                                        &Ywork, &Ework, cm1) ;
                                }
                                t = CPUTIME - t ;
                                t = MAX (t, 0) / NTRIALS ;

                                /* check the solution and log the time */
                                Xsetp = Xset->p ;
                                Xseti = Xset->i ;
                                xlen = Xsetp [1] ;
                                X1x = X->x ;
                                X2x = X2->x ;
                                Lnz = L->nz ;

                                /*
                                  printf ("\ni %d xlen %d  (%p %p)\n", i, xlen, X1x, X2x) ;
                                */

                                if (xtype == CHOLMOD_REAL)
                                {
                                        fl = 2 * xlen ;
                                        for (k = 0 ; k < xlen ; k++)
                                        {
                                                j = Xseti [k] ;
                                                fl += 4 * Lnz [j] ;
                                                err = X1x [j] - X2x [j] ;
                                                err = ABS (err) ;
                                                resid [3] = MAX (resid [3], err) ;
                                        }
                                }
                                else
                                {
                                        fl = 16 * xlen ;
                                        for (k = 0 ; k < xlen ; k++)
                                        {
                                                j = Xseti [k] ;
                                                fl += 16 * Lnz [j] ;
                                                err = X1x [2*j  ] - X2x [2*j  ] ;
                                                err = ABS (err) ;
                                                resid [3] = MAX (resid [3], err) ;
                                                err = X1x [2*j+1] - X2x [2*j+1] ;
                                                err = ABS (err) ;
                                                resid [3] = MAX (resid [3], err) ;
                                        }
                                }
                                if (timelog) fprintf (timelog, "%g %g %g %g\n",
                                                      (double) i, (double) xlen, fl, t);

                                /* clear B for the next test */
                                if (xtype == CHOLMOD_REAL)
                                {
                                        B2x [i] = 0 ;
                                }
                                else
                                {
                                        B2x [2*i  ] = 0 ;
                                        B2x [2*i+1] = 0 ;
                                }

                        }

                        if (timelog)
                        {
                                fprintf (timelog, "] ; resid = %g ;\n", resid [3]) ;
                                fprintf (timelog, "lnz = %g ;\n", cm1->lnz) ;
                                fprintf (timelog, "t = %g ;   %% dense solve time\n", ts [2]) ;
                                fclose (timelog) ;
                        }

                        resid [3] = resid [3] / cholmod_norm_dense (X, 1, cm1) ;

                        cholmod_free_dense (&Ywork, cm1) ;
                        cholmod_free_dense (&Ework, cm1) ;
                        cholmod_free_dense (&X2, cm1) ;
                        cholmod_free_dense (&B2, cm1) ;
                        cholmod_free_sparse (&Xset, cm1) ;
                        cholmod_free_sparse (&Bset, cm1) ;
                }

                /* ------------------------------------------------------------------ */
                /* compute the residual */
                /* ------------------------------------------------------------------ */
                cholmod_write_dense (f, X, NULL, cm1) ;
                
                if (method < 3)
                {

                        if (A->stype == 0)
                        {
                                /* (AA'+beta*I)x=b is the linear system that was solved */
                                /* W = A'*X */
                                W = cholmod_allocate_dense (A->ncol, 1, A->ncol, xtype, cm1) ;
                                cholmod_sdmult (A, 2, one, zero, X, W, cm1) ;
                                /* R = B - beta*X */
                                R = cholmod_zeros (n, 1, xtype, cm1) ;
                                Rx = R->x ;
                                Xx = X->x ;
                                if (xtype == CHOLMOD_REAL)
                                {
                                        for (i = 0 ; i < n ; i++)
                                        {
                                                Rx [i]  = Bx [i] - beta [0] * Xx [i] ;
                                        }
                                }
                                else
                                {
                                        /* complex case */
                                        for (i = 0 ; i < n ; i++)
                                        {
                                                Rx [2*i  ] = Bx [2*i  ] - beta [0] * Xx [2*i  ] ;
                                                Rx [2*i+1] = Bx [2*i+1] - beta [0] * Xx [2*i+1] ;
                                        }
                                }
                                /* R = A*W - R */
                                cholmod_sdmult (A, 0, one, minusone, W, R, cm1) ;
                                cholmod_free_dense (&W, cm1) ;
                        }
                        else
                        {
                                /* Ax=b was factorized and solved, R = B-A*X */
                                R = cholmod_copy_dense (B, cm1) ;
                                cholmod_sdmult (A, 0, minusone, one, X, R, cm1) ;
                        }
                        rnorm = cholmod_norm_dense (R, 0, cm1) ;	    /* max abs. entry */
                        xnorm = cholmod_norm_dense (X, 0, cm1) ;	    /* max abs. entry */

                        axbnorm = (anorm * xnorm + bnorm + ((n == 0) ? 1 : 0)) ;
                        resid [method] = rnorm / axbnorm ;
                }
        }

        tot = ta + tf + ts [0] ;

        /* ---------------------------------------------------------------------- */
        /* iterative refinement (real symmetric case only) */
        /* ---------------------------------------------------------------------- */

        resid2 = -1 ;
        if (A->stype != 0 && A->xtype == CHOLMOD_REAL)
        {
                cholmod_dense *R2 ;

                /* R2 = A\(B-A*X) */
                R2 = cholmod_solve (CHOLMOD_A, L, R, cm1) ;
                /* compute X = X + A\(B-A*X) */
                Xx = X->x ;
                Rx = R2->x ;
                for (i = 0 ; i < n ; i++)
                {
                        Xx [i] = Xx [i] + Rx [i] ;
                }
                cholmod_free_dense (&R2, cm1) ;
                cholmod_free_dense (&R, cm1) ;

                /* compute the new residual, R = B-A*X */
                R = cholmod_copy_dense (B, cm1) ;
                cholmod_sdmult (A, 0, minusone, one, X, R, cm1) ;
                rnorm2 = cholmod_norm_dense (R, 0, cm1) ;
                resid2 = rnorm2 / axbnorm ;
        }

        cholmod_free_dense (&R, cm1) ;
        
        /* ---------------------------------------------------------------------- */
        /* print results */
        /* ---------------------------------------------------------------------- */

        anz = cm1->anz ;
        for (i = 0 ; i < CHOLMOD_MAXMETHODS ; i++)
        {
                fl = cm1->method [i].fl ;
                xlnz = cm1->method [i].lnz ;
                cm1->method [i].fl = -1 ;
                cm1->method [i].lnz = -1 ;
                ordering = cm1->method [i].ordering ;
                if (fl >= 0)
                {
                        printf ("Ordering: ") ;
                        if (ordering == CHOLMOD_POSTORDERED) printf ("postordered ") ;
                        if (ordering == CHOLMOD_NATURAL)     printf ("natural ") ;
                        if (ordering == CHOLMOD_GIVEN)	     printf ("user    ") ;
                        if (ordering == CHOLMOD_AMD)	     printf ("AMD     ") ;
                        if (ordering == CHOLMOD_METIS)	     printf ("METIS   ") ;
                        if (ordering == CHOLMOD_NESDIS)      printf ("NESDIS  ") ;
                        if (xlnz > 0)
                        {
                                printf ("fl/lnz %10.1f", fl / xlnz) ;
                        }
                        if (anz > 0)
                        {
                                printf ("  lnz/anz %10.1f", xlnz / anz) ;
                        }
                        printf ("\n") ;
                }
        }
        
        printf ("ints in L: %15.0f, doubles in L: %15.0f\n",
                (double) isize, (double) xsize) ;
        printf ("factor flops %g nnz(L) %15.0f (w/no amalgamation)\n",
                cm1->fl, cm1->lnz) ;
        if (A->stype == 0)
        {
                printf ("nnz(A):    %15.0f\n", cm1->anz) ;
        }
        else
        {
                printf ("nnz(A*A'): %15.0f\n", cm1->anz) ;
        }
        if (cm1->lnz > 0)
        {
                printf ("flops / nnz(L):  %8.1f\n", cm1->fl / cm1->lnz) ;
        }
        if (anz > 0)
        {
                printf ("nnz(L) / nnz(A): %8.1f\n", cm1->lnz / cm1->anz) ;
        }
        printf ("analyze cputime:  %12.4f\n", ta) ;
        printf ("factor  cputime:   %12.4f mflop: %8.1f\n", tf,
                (tf == 0) ? 0 : (1e-6*cm1->fl / tf)) ;
        printf ("solve   cputime:   %12.4f mflop: %8.1f\n", ts [0],
                (ts [0] == 0) ? 0 : (1e-6*4*cm1->lnz / ts [0])) ;
        printf ("overall cputime:   %12.4f mflop: %8.1f\n", 
                tot, (tot == 0) ? 0 : (1e-6 * (cm1->fl + 4 * cm1->lnz) / tot)) ;
        printf ("solve   cputime:   %12.4f mflop: %8.1f (%d trials)\n", ts [1],
                (ts [1] == 0) ? 0 : (1e-6*4*cm1->lnz / ts [1]), NTRIALS) ;
        printf ("solve2  cputime:   %12.4f mflop: %8.1f (%d trials)\n", ts [2],
                (ts [2] == 0) ? 0 : (1e-6*4*cm1->lnz / ts [2]), NTRIALS) ;
        printf ("peak memory usage: %12.0f (MB)\n",
                (double) (cm1->memory_usage) / 1048576.) ;
        printf ("residual (|Ax-b|/(|A||x|+|b|)): ") ;
        for (method = 0 ; method <= 3 ; method++)
        {
                printf ("%8.2e ", resid [method]) ;
        }
        printf ("\n") ;
        if (resid2 >= 0)
        {
                printf ("residual %8.1e (|Ax-b|/(|A||x|+|b|))"
                        " after iterative refinement\n", resid2) ;
        }

        printf ("rcond    %8.1e\n\n", rcond) ;

        if (L_is_super)
        {
                cholmod_gpu_stats (cm1) ;
        }
        /*
        FILE* f3;
        if ((f3 = fopen (argv[3], "r")) == NULL)
        {
                my_handler (CHOLMOD_INVALID, __FILE__, __LINE__,
                            "unable to open node index file") ;
        }
        
        char nodes[n][20];
        
        fgets(line, 65532, f3);
        temp = strtok(line, "\t\n");
        strcpy(nodes[0], temp);
        //strncpy(nodes[0], temp, strlen(temp)+1);
        for(i = 1; i < n; i++){
                printf( "h %d\n", i);
                temp = strtok(NULL, "\t\n");
                //strncpy(nodes[i], temp, strlen(temp)+1);
                strcpy(nodes[i], temp);
        }
        */
        cholmod_write_dense (f__, B, NULL, cm1) ;
        printf("X size is %d %d\n", X->nrow, X->ncol);
        cholmod_write_dense (f, X, NULL, cm1) ;

        cholmod_free_factor (&L, cm1) ;
        cholmod_free_dense (&X, cm1) ;
        
        /* ---------------------------------------------------------------------- */
        /* free matrices and finish CHOLMOD */
        /* ---------------------------------------------------------------------- */

        cholmod_free_sparse (&A, cm1) ;
        cholmod_free_dense (&B, cm1) ;
        cholmod_finish (cm1) ;
    
        return (0) ;
}
